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Theoretical studies of the production of real thermal photons in relativistic heavy ion collisions at 
the Relativistic Heavy Ion Collider (RHIC) are performed. The space-time evolution of the colliding 
system is modelled using MUSIC, a 3+1D relativistic hydrodynamic simulation, using both its ideal 
and viscous versions. The inclusive spectrum and its azimuthal angular anisotropy are studied 
separately, and the relative contributions of the different photon sources are highlighted. It is shown 
that the photon V2 coefficient is especially sensitive to the details of the microscopic dynamics like 
the equation of state, the ratio of shear viscosity over entropy density, r//s, and to the morphology 
of the initial state. 



The study of relativistic collisions of nuclei constitutes a vibrant branch of subatomic physics that straddles nuclear 
and particle physics. It offers a privileged window on the physics of hot and dense strongly interacting matter and 
as such, it complements astrophysical studies. There, the hadronic equation of state is an ingredient of paramount 
importance that enters the evaluation of the bulk properties of neutron stars, for example. In comparison, relativistic 
nuclear collisions do offer the considerable practical advantage of providing laboratory control over the projectile and 
target characteristics, together with the beam energy. This physics currently defines a large experimental effort being 
pursued at the Relativistic Heavy Ion Collider (RHIC), at Brookhaven National Laboratory and, more recently, at 
the Large Hadron Collider (LHC), at CERN. One of the remarkable results that emerged from the RHIC program 
so far, is the fact that the hot and dense hadronic matter produced there [HI] could be described using almost ideal 
hydrodynamics, that is with a small shear viscosity coefficient rj [5J, compared to the entropy density s. The first LHC 
flow results for heavy ion collisions [5] also suggest similar conclusions: a recent overview of flow results can be found 
in Ref. [7]. In fact, the progress in both theoretical and in experimental analyses has been such that the goal of a 
quantitative extraction of the shear viscosity coefficient of hot and dense strongly interacting matter from relativistic 
nuclear collision data now appears closer than ever [8lU3|. 

In heavy-ion collisions, the flow has been characterized by considering a Fourier expansion of the the triple differential 
cross-section, with the variable being the azimuthal angle with respect to the reaction plane [14] : 



where ip r is the reaction plane angle. The expansion coefficients, v n , will then quantify the degree of azimuthal 
anisotropy. In the progress towards a precise and quantitative characterization of the hydrodynamical state of nuclear 
collisions, a recent development consisted of linking the odd-numbered coefficients to fluctuations in the initial state 
[15] . Indeed, up to this point the flow observables in nuclear collisions had been analyzed by considering the coefficients 
of the Fourier expansion of the azimuthal angle distribution of the particle spectra, together with smooth initial 
conditions for which the odd values of the Fourier coefficients vanish [14] . It is fair to write that a rich and quantitative 
picture of nuclear flow is now emerging] 

In general, the measured hadronic observables give a dynamical snapshot of the conditions that existed on the 
freeze-out hypersurface. In contrast, electromagnetic radiation is emitted throughout the space-time evolution and 
suffers negligible final-state interactions, owing mainly to the smallness of a, the electromagnetic coupling constant. 
Real and virtual photons are thus penetrating probes, and as such can carry information about the different stages 
of the high energy collisions. A consequence of this statement is that accurate and meaningful calculations of photon 
spectra in relativistic nuclear collisions will need realistic electromagnetic emissivities and precise modelling of the 
space-time dynamics. The goal in this paper is to extend the calculations of real photon production to situations 
which incorporate the developments made on the purely hadronic front. Cases where the emitting source is no longer 
in local thermal equilibrium will be considered, together with cases where the initial states of the nucleus-nucleus 



A discussion on how to calculate the elliptic flow, V2, in the presence of fluctuating initial conditions appears later in this paper. 
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collisions are no longer smooth but are allowed to fluctuate event-by-event. The paper is organized as follows: section 
[il] contains a brief description of MUSIC, our implementation of 3+1 viscous hydrodynamics which is used to calculate 
the evolution of the background medium. In section [nil we gi ye a short explanation of viscous corrections to the local 
momentum distribution function and the photon emission rates used in this study. Our main results are presented in 
section |IV] and we conclude in section |Vl 

II. HYDRO DYNAMICAL EVOLUTION 

As mentioned earlier, photons are penetrating probes that are emitted throughout the heavy-ion collision. It 
is thus imperative to evaluate their observed properties with a time-evolution scenario that is both realistic and 
consistent with a large number of empirical observables. One such approach is MUSIC, a three-dimensional simulation 
of relativistic hydrodynamic systems [16]. The general features of MUSIC are described below, first in the ideal limit , 
and then incorporating a finite coefficient of shear viscosity. 

The solution of the conservation laws for the stress-energy tensor and the net baryon current, T^ v and Jg, respec- 
tively, dictate the evolution in time of an ideal hydrodynamical system. More specifically, 

«L = o, 

«,idcal = (2) 

and 

4 l ,idcal=PB^ (3) 

Note that P is the local pressure, e is the local energy density, pb is the local net baryon density, u M = (7, 7V) is the 
local flow velocity with respect to some fixed frame, and g^ v = diag (1,-1,-1,-1). This represents a set of five scalar 
equations, with six unknowns. The set is closed by specifying an equation of state, P(£,/?b)- MUSIC is implemented 
in r — r\ s coordinates, where r is the proper time, and r\ s the space-time rapidity. The transformations to real time 
and longitudinal coordinate variables, {£, z}, are 

t = t cosh rj s , z — t sinh r] s (4) 

The solution of Eqs. ^ in t — rj s coordinates is obtained with the Kurganov-Tadmor method [17] , and an equation 
of state extracted from lattice QCD calculations [H] is used in this work. Importantly, MUSIC is a three-dimensional 
simulation, and is therefore capable of following the time evolution of the rapidity profile. It has been used for the 
successful calculation of flow variables, including elliptic flow and higher flow harmonics [16] . A discussion of results 
with ideal hydrodynamics is postponed, turning now to the inclusion of viscous effects. 

The first-order - or Navier-Stokes - formalism for viscous hydrodynamics is known to introduce unphysical super- 
luminal signals that spoil the theory's stability. Various formulations of second-order hydrodynamics [19-23 address 
this problem, and a variant |24j of the Israel-Stewart formalism is used here. In this approach, the stress-energy tensor 
is = T(^ al + tt^, and the evolution equation are 

A£A£u CT <9 CT 7r Q/3 = - — (tt^ - S^) - %tx^ (d a u a ) 

Trr 3 

(5) 

where A MI/ = g^ v — u^u v , and tv is usually interpreted as a relaxation time . The first-order (in velocity gradients) 
viscous part of the stress-energy tensor appears here and is 

= 7? ^A"u v + A v u^ - ^""Vo^j , (6) 

with the coefficient of shear viscosity 77, and V M = A^ u d u . The viscous stress-energy tensor tt^ u is clearly a complicated 
object that is evaluated dynamically. Finally, vorticity and numerically small terms have been neglected. 

Hydrodynamic calculations require their initial states to be defined. In this work, smooth or averaged initial 
conditions (AIC) will be considered, as well as cases where these initial states are allowed to fluctuate (FIC) about 
that average. The procedure for implementing AICs in MUSIC is described in detail in Ref. [TCj, and that for FICs, 
in Ref. [5]; those descriptions will not be repeated here. The initial time for the hydro to start is defined by tq, in 
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this work this value is set to To = 0.2 fm/c, and the freeze-out energy density is 0.12 GeV/fm 3 , which approximately 
corresponds to T = 137 MeV. 

Considering first AICs, it is instructive to study how the inclusion of shear viscosity affects the bulk evolution. The 
physical case being considered is that of Au + Au, at y/s — 200 GeV, at an impact parameter of b = 4.47 fm, which 
represents a - 20 % centrality class. As mentioned earlier, the equation of state used here is the parametrization 
"s95p-vl" from Ref. [TH]. In this parametrization, the fit to the lattice QCD data is made above T = 250 MeV with 
the constraint that at T — 800 MeV, the energy density reaches 95 % of the Stefan-Boltzmann value. It is also worth 
noting that this parameterization correctly reproduces the trace anomaly around the transition temperature: see Ref. 
[18] for more details. 

Figure [T] shows the evolution of temperature for a fixed cell at x — y = 2.5 fm and z = 0. In the case with a non-zero 
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FIG. 1. (Color online) Left panel: The evolution of temperature in a fixed cell, as a function of time. Results for ideal and 
viscous hydrodynamics (with r]/s = l/47r) are shown. Right panel: A closer view of the early time temperature evolution. 



shear viscosity, a value of r\js (the shear viscosity divided by the entropy density) = 1/Aix has been used. This value 
has been suggested as a lower universal bound [25], a statement which has raised some controversy and thus needed 
to be qualified [2T>H2"S] . Furthermore, r]/s will also depend on the local temperature of the medium, but a constant 
value will suffice for the study in this work. Note that the finite viscosity calculation has a smaller initial temperature 
than the ideal one, as entropy will grow in the viscous case, affecting the observed final particle multiplicity: it is 
important to compare calculations consistent with a given set of hadronic observables. 



III. PHOTON EMISSION FROM IDEAL AND VISCOUS MEDIA 

Viscous corrections on microscopic processes involving particles have been included by writing the in-medium 
distribution functions with an out-of-equilibrium correction, / — > /q + Sf, where / is an ideal Bose-Einstein/Fermi- 
Dirac distribution function. This is most easily seen by considering the particle spectra being generated from the 
Cooper-Frye formalism [29], and requiring that the energy momentum tensor be continuous across the freeze-out 
hypersurface. In a multispecies ensemble, a popular ansatz that satisfies the continuity requirements is 

6 ft = foi (1 ± /oi^V ^* (7) 

for the distribution function of species i. This form is used in this work. In general, there can be an overall constant, 
different for each species, that multiplies Eq. ^ j2U]. It is implicit in this treatment that Sf should represent a small 
correction. In the calculation of pions and of other hadronic observables, one simply needs to verify this statement on 
the freeze-out hypersurface. This is not the case for electromagnetic emission which occurs at all time-scales of the 
hydro evolution. Therefore, photon calculations in viscous media will represent a stringent test of the validity of the 
viscous dynamics, as shall be seen later. 
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FIG. 2. The Compton and quark- ant iquark annihilation contributions to photon production. 



A. Photon emission from the QGP 

Rates, complete at leading order in a s , for the emission of photons from a thermal ensemble of partons have 
now been available for a decade [3T]. The extension of these results to viscous media necessitates revisiting the 
resummation procedure in Ref. 31j with out-of-equilibrium distributions: a process we shall not perform here. We 
rather concentrate on a subset of the diagrams: the Compton and quark-antiquark annihilation processes shown 
in Figure [2j It is instructive to compare the photon rate obtained through the approach described above with the 
complete result at leading order in a s : this is done in Fig. [3| At low px, the full leading order rates are an order of 
magnitude larger than the naive leading order rates owing to additional processes. For examples, the former receive 
a large contribution of brehmsstrahlung from quarks of all momenta in this range. For pr > 1 GeV, the full leading 
order rates are only larger by about a factor of two (there is however some temperature dependence to the position 
of this transition window). 
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FIG. 3. (Color online) A comparison of the equilibrium photon rate from the processes shown in Figure [5] (dashed lines) with 
that obtained tallying all channels contributing at leading order in a s (full lines), for Nf — 3. The lower set of curves are for 
T = 250 MeV, and the upper ones are for T = 350 MeV. 



The net photon emission rate R, summing these individual processes of the type 1 + 2 — > 3 + 7, is obtained by 
evaluating 
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where the coefficients a,b,c are defined in Eq. (A. 11), and where |A^i| 2 = I6irs 2 d<7i/dt, with 

87raa s u 2 + t 2 
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Note also the degeneracy factors jVannihil. = 20, and Acompt. = 320/3, for Nf = 2. In the case of Nf — 3, those 
numbers become A/" a nnihil. = 24, and Acompt. = 384/3. Calculations of the photon production rate from these channels 
were done in Ref. [32j . an evaluation with general anisotropic distribution functions (not limited to small deviations 
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from equilibrium) appeared in Ref. |33] , and a viscosity-corrected rate (to first order in df) was obtained recently in 
[34j , assuming forward-scattering dominance of the photon-producing reaction. The rates reported here are obtained 
through a numerical integration of Eq. Q with out-of-equilibrium distribution functions (Eq. |7])). The integrations 
span the entire accessible phase space, carefully avoiding divergences as prescribed in Ref. [32 . Appropriate quantum 
statistics have been used. 



B. Photon emission from the hadronic gas 



As the ensemble of partons thermalizes (totally or partially) and then expands and cools, it hadronizes into an 
ensemble of colorless hadrons called here the hadronic gas (HG) which continues to expand and to cool even more. 
The HG thermal electromagnetic emissivity has been characterized in Ref. [35 . Following that reference, a Massive 
Yang-Mills (MYM) model is used to model the interactions between light pseudoscalars, vector and axial vector 
mesons. The set we consider contains the elements {77, K, p, K* , 01}, and the most important photon-producing rates 
are 7r + p — s- 7T + 7, tt + tt — > p + 7, ir + K* — > K + 7, n + K — > K* + 7, p + K —> K + 7, K* + K — >• tt + 7. 
Two-body photon-production processes dominate the phase space for photon transverse momenta above 0.5 GeV [35] . 
All isospin-allowed channels are considered. 

The viscous corrections also demand a complete recalculation of the HG photon rates, by including the corrected 
distribution functions - see Eq. Q - in all the relevant rate equations. Note that corrections of order 5f 2 are neglected 
for consistency, as are corrections to Pauli-blocking or Bose-enhancement effects. These corrections are found to be 
small. The Appendix outlines the procedure for correcting the electromagnetic emissivities, allowing for viscous effects 
in the hadronic distribution functions. 



IV. RESULTS 



A. Viscous corrections: generalities 



For both cases discussed in the previous section (QGP and HG), rates for "viscous photons" were not shown. In 
fact, those require detailed dynamical information as they depend on the details of tt' 11 ' and of its time evolution as 
specified by Eqs. ^ and ([5]). It is thus appropriate to examine this quantity here, and this is done in Figure [4j in 
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FIG. 4. (Color online) Left panel: The time evolution of different components of the local -e^ v tensor, divided by r\. Right 
panel: The time evolution of the diagonal elements of (scaled by -q), and also that of the trace of the viscous tensor. The 
calculations are done for a fluid cell at x — y = 2.5 fm, and 2 = 0, and the impact parameter is b = 4.47 fm. 

the rest frame of a fluid cell; note that there tt u is 0. At the initial time, the viscous corrections are non-existent, 
as we initialize the viscous pressure tensor to zero. They build up quickly, and then decay back to zero. Right after 
the initial time, the magnitude of the zz component is larger than the other two diagonal ones by roughly a factor of 
2, and this fact persists up to late times. The relative sign of tt zz can be understood from the fact that "k^ should 
be traceless in the fluid rest frame (c.f. Eqs. (|5]|6|). Note that this requirement was not enforced explicitly at each 
step of the calculation. The preservation of this trace then reflects the stability of the numerics: see the right panel 



() 



of Figure |4j The slight difference between ir xx and Tr yy is to be expected because of the elliptic shape of the system: 
the x — y symmetry is broken by the finite impact parameter. 
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FIG. 5. (Color online) The fraction of rj s ~ fluid cells with a certain value of Sf / fo, for different values of the photon momentum 
in the nucleus-nucleus center of mass frame: 2 GeV (left panel), and 3 GeV (right panel). The range of temperature T > 250 
MeV corresponds to t — To < 2 fm/c. 



To get a qualitative picture and develop some intuition for the importance of the viscous corrections, the following 
procedure was implemented. In the nucleus-nucleus center of mass frame, one picks a photon momentum in the x — y 
plane, at an angle of 7r/4. The z axis is the beam axis. Lorentz-transforming to the rest-frame of each fluid cell, 
the local value of the photon momentum is obtained. Since photons are formed in 2 — > 2 processes, the magnitude 
of this momentum is then roughly equal to the magnitude of the momentum of one of the interacting particles. 
Finally, posing that this particle is a massless fermion will enable a determination of the viscous correction to its 
distribution function. This study is restricted to a slice in space-time rapidity, r] s , centered around 0. This procedure 
is clearly approximate - not all particles are massless and not all are fermions - but should nevertheless produce a 
result indicative of the physics at play. Defining bins of size 0.1 in the relative variable Sf / fo, the fraction of j] s w 
cells with a certain value of this relative variable is plotted in Figure [5] In addition, in this study we concentrate on 
two ranges of temperature: one corresponding to "early times" , and another corresponding to "late times" according 
to Figure [T] The photon energies chosen are typical values of the photon spectrum, see the next subsection. For 
a photon energy of 2 GeV, one sees that ~ 20% of the fluid cells have a Sf/fo > 1, at early times, and that the 
distribution around this value is fairly narrow. For a higher photon energy of 3 GeV, this distribution has grown in 
width, now with 80% of the high-temperature cells with Sf/fo > 1 and ss 30% of them with Sf/fo > 2: a clear 
violation of the perturbative nature of the approximation. In a given panel of Figure [5j the amount of larger viscous 
corrections at higher temperatures can be understood: those cases correspond to situations at early times where the 
elements of the it^ u tensor are large (see Figure Q). Higher momenta will command larger viscous corrections, see Eq. 
Q , and hence the broadening of the distributions for a given range in T when going from the left to the right panel of 
Figure [5j Note that repeating the lower temperature part of this analysis by assuming that the corrected particles are 
massive 7r's or even massive p's do not change its conclusions. A negligible but non-zero number of high-momentum 
cells have Sf/fo < — 1- For those cells the photon emission probability has been set to zero. What emerges here is 
an explicit message of caution. In contrast to hadron calculations where only properties at the freeze-out surface are 
required, the evaluation of electromagnetic signals requires the entire time evolution to be monitored. In the case of 
the analysis shown here, one should also keep in mind that thermal photons with an energy equal to, or greater than, 
3 GeV will lie below those from other sources like the direct photons from primordial nucleon-nucleon collisions, for 
example [551 157] . 

One now proceeds to the evaluation of photon characteristics, and of the influence of viscous effects on them. All 
of the calculation results shown in the next sections will rely on the use of Averaged Initial Conditions (AICs). A 
discussion of the effects of Fluctuating Initial Conditions (FICs) on real photons first appears in section IV D 



B. Photon spectrum 



As mentioned previously, the measured photon spectrum receives contributions from all times and all phases spanned 
by the collision dynamics. There are treated in turn here, for pedagogical purposes. The photons originating solely 
from the QGP phase are first shown in Figure [6j and they are obtained by integrating the rate in Eqs. ^ and (A.26) 
throughout the time evolution dictated by MUSIC. 
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FIG. 6. (Color online) Left panel: The photon yield originating from the phase with parton degrees of freedom only. The 
contribution from ideal hydro is shown (solid curve), together with the result of using a time evolution associated with viscous 
hydrodynamics (dotted line), and using a viscous time evolution and corrected microscopic distribution functions (dash-dotted 
line). Right panel: The photon yield originating from the hadronic gas only. The meaning of the different curves is the same 
as that in the left panel. 



As compared to an ideal hydrodynamical evolution, the viscous evolution starts with a lower initial temperature 
when the system is entirely in its QGP phase (see Figure [lj. Therefore, integrating the QGP photon rates with 
a viscous hydrodynamic evolution alone produces a photon spectrum slightly lower at high values of px than that 
generated by the ideal hydro. This is the dotted line in Figure [6] Then, using the corrected distribution functions 
make the photon spectrum harder, as the correction grows as a function of momentum. This is the upper curve in 
the same figure. The hardening of the QGP photon spectra owing to shear viscous effects had also been noticed in 
previous work [34 ) [38 l 139] . 

Turning now to photons originating solely from the HG sector, the relevant spectrum is shown in the right panel of 
Figure [6] Interestingly, the spectrum with the viscous corrections (viscous hydro and corrected distribution functions) 
is essentially undistinguishable from that obtained using ideal rates integrated with a viscous time evolution. This can 
be understood by considering the fact that photons from the HG are emitted later in time, essentially when 7r M " ~ 0, 
as is made clear in Figure[4j The effect of viscosity are manifested in a slightly harder spectrum: in part a consequence 
of the temperature in the viscous evolution remaining higher than that in the ideal evolution for intermediate and late 
times, as shown in Figure[l] The yield of real photons from all thermal sources (QGP + HG) is shown in FigureJTJ for 
an ideal hydrodynamic evolution and also for a viscous evolution (viscous hydro and corrected distribution functions) . 
The difference between the two scenarios is actually small at intermediate values of the photon transverse momentum, 
growing to being approximately 100% at px = 4 GeV. At that energy however, the purely thermal photons will lie 
below other sources and will be sub-dominant, as mentioned already. One may thus conclude here that extracting 
information about the shear viscosity from photon spectra alone will be an arduous task. More work is needed however 
to include all the photon sources in a theoretically consistent way, with all the viscous corrections. 



C. Photon elliptic flow 



The flow characteristics of hadrons have contributed considerably to quantify the details of the underlying hydro- 
dynamics, and this fact has been hailed as one of the major milestones of the RHIC program. As for photons, their 
elliptic flow holds the potential of providing more insight into the dynamics of heavy ion collisions, and into the phase 
structure of QCD. Indeed, the shape of the real photon w 2 coefficient is directly sensitive to the nature of the under- 
lying degrees of freedom jJU], unlike the single-photon spectra of the previous section. The elliptic flow of photons 
originating solely from the QGP is shown in Figure [8] If one neglects the correction to the distribution functions, 
the elliptic flow from the viscous evolution appears slightly larger than in the ideal case, reflecting an increase in 
the azimuthal asymmetry of the fluid flow pattern due to viscosity, consistent with the large gradients at early times 
implied by Figure [4] However, the corrections to the distribution functions dominate, and make the net anisotropy 
even smaller than in the ideal case. This behaviour is consistent with the results of Ref. [39]. It is also worthwhile to 
point out that the apparently negative values of %}% at very low momenta for photons in the QGP are obtained only 
using the photon-production rates corresponding to Figure [2] The complete rates of Ref. [31] do not yield negative 
photon elliptic flow in ideal hydrodynamics. The negative values also appear in earlier calculations [311 [Ml- The HG 
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FIG. 7. (Color online) The net thermal photon yield, from QGP and HG sources. The ideal spectrum (i.e. using an ideal 
hydrodynamics background), and the viscous spectrum (using a viscous hydrodynamics background and corrected microscopic 
distribution functions) are shown as a solid and dotted line, respectively. 
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FIG. 8. (Color online) Left panel: The thermal photon elliptic flow, considering only the photons originating from the QGP. As 
in previous figures, the results of using ideal hydrodynamics (solid line), viscous hydrodynamics with equilibrium rates (dotted 
line), and viscous hydrodynamics with Sf corrections (dash-dotted line) are shown separately. Right panel: The thermal photon 
elliptic flow, considering only the photons originating from the HG. The lines have the same meaning as those in the left panel. 

V2 is shown in the right panel of Figure [8] and there, all viscous corrections make the elliptic flow smaller, unlike the 
case for the QGP. This is again a reflection of the richness of the dynamics contained in the time-dependence of ir^". 
Further note that the small structure at low momenta signals a crossover between two different hadronic channels 
[4"tl] . The net photon w 2 is then calculated and shown in Figure [9j Importantly, the total v-i is a weighted average of 
the individual (QGP, and HG) coefficients, the weight being the value of the appropriate single-photon distribution. 
Hence, in the computation of the final V2, the small QGP v-i will get multiplied by a large emission rate, whereas 
the smaller emission rate of the HG phase gets partially compensated by the larger flows. Both phases therefore 
contribute to the final profiles shown in Figure [9] 

D. Fluctuating initial conditions (FIC) 

The recent years have witnessed a paradigm-shift in the analysis of heavy ion collision data. Up until recently, 
smooth initial state distributions were mostly used in hydrodynamics analyses of relativistic nuclear collisions. These, 
together with conservation laws, imply that odd-numbered expansion coefficients in Eq. ^ vanish identically. As 
discussed in the Introduction, this situation has changed with the work of Ref. [15] linking odd-numbered flow 
harmonics to initial state fluctuations. The hydrodynamic simulation MUSIC with viscous corrections has recently 
been modified to include FICs [S]. This has been used to make a prediction for size and momentum dependence of 
the hadronic U3 at RHIC. This prediction has been recently confirmed 0T]. Here we seek to assess the importance of 
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FIG. 9. (Color online) The net thermal photon elliptic flow. The curves have the same meaning as in Figure [7] 

the event-by-event fluctuations on photon observables. 

For initial conditions that are not smooth, it is important to specify how the reaction plane is determined. The 
"participant plane" is used here. Namely, one calculates event-by-event the angle ip2 with respect to the reaction 
plane defined by the impact parameter: 
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(10) 



where the averages are over wounded nucleon positions, {r,4>), in the transverse plane. The angle ip2 then goes into 
the evaluation of u 2 , with ip2 replacing if> r in Eq. ([I]). Note that the initial eccentricity is maximized by the choice of 
this participant plane. The studies performed here used ensembles of 50 events, leading to uncertainties of the order 
of 5% on thermal photon spectra, and of the order of 15% on thermal photon v^- The precise value of these variations 
is of course p^-dependent, but we find that elliptic flow does depend more strongly on the initial structure of the 
energy density distribution than the momentum spectrum. 

As already observed for hadrons [H] and more recently for photons |H], the lumpy initial states lead to a yield 
enhancement. Again, the QGP and HG contributions are calculated separately. They are shown in the two panels 
of Figure |10[ and the quantitative importance of the enhancement can be judged there. As done previously, only 
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FIG. 10. (Color online) The thermal photon yield, showing the effect of FICs. The left panel shows the contribution from the 
QGP, the right panel that of the HG. Note that the curve labeled "FIC" also includes all viscous corrections (time evolution 
and 5f) 

this time with FICs, we plot the thermal photon i>2 for QGP and HG. This is shown in Figure [TT] Finally, the net 
photon spectrum and V2 are shown in Figure |12| Clearly, in the centrality range studies in this work, the hot spots 
and large gradients generated by the fluctuating initial conditions lead to a harder photon spectrum and to a larger 
elliptic flow, and this remains true with the inclusion of a finite shear viscosity to entropy density ratio. 
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FIG. 11. (Color online) The thermal photon V2, showing the effect of FICs. The left panel shows the contribution from the 
QGP, and the right panel that of the HG. Note that the curve labeled "FIC" also includes all viscous corrections (time evolution 
and 5f). 
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FIG. 12. (Color online) The net thermal photon yield (left panel) and «2 (right panel), showing the effect of FICs. Note that 
the curve labeled "FIC" also includes all viscous corrections (time evolution and 5f). 



V. CONCLUSION 



In this work we have sought to establish the quantitative importance of a finite shear viscosity coefficient and of 
fluctuating initial conditions on two real photon observables: the one-body spectrum and the transverse momentum 
dependence of the elliptic flow coefficient. This was done using MUSIC, a realistic 3+1D relativistic hydrodynamical 
simulation. Importantly, comparisons between cases with and without viscous corrections were done using conditions 
tuned to hadronic experimental data, and this was the case also for studies involving FICs. Results obtained here show 
that the combined effects of the viscosity and of the FICs are large enough to make their inclusion mandatory in any 
attempt to quantitatively extract transport coefficients of the hot and dense matter from thermal photon data. It was 
not the point of this work to explicitly compare with experimental measurements just yet. Firstly, 3+1D relativistic 
viscous hydrodynamics models are in their infancy, and systematic studies of all parameter dependences, in the spirit 
of that in Ref. [15] for example, will be useful to establish a more precise quantitative link between observables and the 
underlying hydrodynamics. Secondly, in what concerns the photon sources, an inclusive and consistent treatment of 
all of them (pQCD photons, photons from jets interacting and fragmenting while losing energy . . . ) with and without 
viscosity is still to be done. Finally, exploring the consequences of what has been found here on electromagnetic 
observables at the LHC should prove interesting and relevant. 

In closing, it is worth mentioning that recently the PHENIX collaboration at RHIC has extracted a direct photon 
V2 from measured data |46) . Interestingly, this analysis concludes that the direct photon elliptic flow is comparable 
in magnitude to that of the 7r°. This large photon elliptic flow is a challenge to most approaches, but may contain 
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some clues about early dynamics prior to hadronic freeze-out [37] ■ However, a complete theoretical analysis including 
all that is known in the literature about the production of direct photons in relativistic heavy ion collisions (viscous 
effects, FICs, the effect of realistic 3+1D hydrodynamical modelling, hadronic chemical potentials, primordial flow) 
is needed before a more precise assessment can be obtained. 
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Appendix: Viscous photon rates 

Our starting point for the photon emission rate of a (1 + 2^3 + 7) process is the expression leading to Eq. Q: 



E 



d 3 R 



d 3 p\ d 3 p 2 d 3 p 3 2n 



d 3 p J 2E 1 (2ir) 3 2E 2 (2it) 3 2E 3 (2tt) 3 2 

x 6\pi +P2-P3- p)/(pi)/(pa) (1 + fips)) 



(A.l) 



The subscripts 1,2 and 3 refer to the two incoming particles and the outgoing particle (not the photon) respectively. 
Notice that the momentum distribution functions depend on the four-momentum and not only on the zero component 
(the energy). This can be written as (with s = (pi + p 2 ) 2 , t = (p 2 — p) 2 , u = —s — t + m 2 + m 2 + m 2 ) 



E 



d 3 R 



dt / dEi / dE 2 / d<f>i / dfc\M\ 



d 3 p 25(2^£2 Jsmin Jtmm jE?ixl J£ „ in JQ Jo 

S( s -m\- m\ - 2E 1 E 2 + 2\ Pl \\p 2 \ cos0)/(px)/(p 2 )(l + f( Pl +p 2 -p)) 



(A.2) 



where 



cos t) = cos oi cos t) 2 + sin tq sin v 2 cos {<p 2 — <p\) 

-t-s + ml + 77il + 2EE 1 „ -t-m% + 2EE 2 

COS6*! = \ , - , 3 , 008^2 = - 

1E\p x \ 2E\p 2 \ 
The 6 angles represent the angle between a particle's three momentum and the photon's. In other words, 

Pa -P 7 = \Pa I \fi<y I COS a 



(A.3) 



(A.4) 



The <p angles represent the azimuthal angles of the incoming particle's direction around the photon's direction. The 
integration boundaries are given by 



s min > (mi +m 2 ) 2 , s min >ml 



(A.5) 



f nin(max) = ^ + ^ 



-b+s/l? 



-( + ) 



1 / s + m\ — m\ \ / s + m 



(s + mf — m 2 ) 2 — 4smf f s — m 



Ef 



2^ 

Em\ m\ — u 
m\ — u 4E 



<E 2 < 



-b-s/b 2 



ac 



(a is negative) 



(A.6) 
(A.7) 
(A.8) 
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The a, b and c coefficients are given below. Now, we can take care of the fa integration with the Dirac delta. Its 
roots are 

_ 1 f-s + mf + m 2 + 2E X E 2 - 2|p 1 ||p 2 1 cos 9 t cos 9 2 \ 

fa=fa± cos . . (A.9) 

V 2|pi||p 2 | sm£>i sin0 2 ' 

However, now we need to keep in mind that where we had f(E 2 ) for the ideal (no viscous corrections) case, we now 
have f(E 2 , 62, fa), so the fa integration goes like 



\ dfaS(s - m\ - m\ - 2E 1 E 2 + 2\ Pl \\p 2 \ cos9)f(E 2 ,9 2 ,fa) 
Jo 

5(fa-fa)f(E 2 ,9 2 ,fa) 
2|pi||p 2 | sin 0i sin# 2 sin (fa - fa) 



i ( — 

sin 0i sin 9 2 \ ./\ 



f{E 2 ,9 2 ,fa) , f{E 2 ,9 2 ,fa) 



2|pi||p 2 | sin 6»i sin6» 2 V \ - C os 2 (0+ - fa) ^1 - cos 2 (0_ - fa) 

■(f(E 2 ,9 2 ,fa) + f(E 2 ,9 2 ,fa) 



2\pi | |p 2 1 sin 9i sin 6> 2 - cos 2 (</>+ - fa) 

E (f(E 2 ,02, fa) + f(E 2 , 9 2 ,fa-)) (A.IO) 



^a~E^+2~b~E 2 
where the a, b and c coefficients are given by 



a = — (s + 1 — m 2 — m 2 ) 2 

b = E((s +t-m\- m|)(s - ra\ - m - 2 2 ) - 2m\(m\ - t)) 

+Ex{ml -t)(s + t-m\ - m 2 ) 
c = c 2 E\ + ciEi + c 
c 2 = -(t - m\) 2 

c\ = — 2E(2m 2 (s + t — m\ ~ m 2 ) — (m 2 — t)(s — m 2 — m 2 )) 

A rp2 2 2 2/ . , 2 3 \ 2/2 ,\2 tt.2 / 2 2\2 

Co = 4_fc/ m 1 m 2 + m 2 (s + t — m 2 — m 2 ) + m 1 (m 2 — t) — E (s — m 1 — m 2 ) 

+(s -m\ - m 2 2 )(t - m 2 2 )(s + t-m 2 2 - m 2 ) (A. 11) 

So at this point, the rate is given by 

(PR 1 C°° /-* mi * x poo /-BJ"™ p2ir 

E ^=w^L ds L dt L?» dEi L» dE2 l ^i 2 /^>^) 

(f(E 2 ,9 2l fa) + f(E 2 ,9 2 , fa)) (1 + f(E 1 +E 2 - E))^==L= = (A.12) 
V / yJaE 2 + 2bE 2 + c 

We have not corrected the distribution functions for the outgoing particle (Pauli blocking or Bose enhancement). At 
this point, we look more carefully at the form of the correction (choosing Bose enhancement for illustrative purposes): 



/(p a ) = f (E a ) + Sf(p a ) = f (E a ) + (^/o(£a)(l + fo(E a ))p a a p^) (A.13) 



One needs to properly express p a in terms of integration variables. Because the incoming particles angles are given 
with respect to the photon's direction, the angles of the photon (in the fluid frame) will explicitly appear in the 
expressions for the particles four momentum. The zero component of p a is just the energy, so the interesting parts 
are the 1,2 and 3 components. In what follows, the subscript a will refer to the incoming particles (1 or 2), and the 
subscript 7 will denote angles for the photon. The photon angles are given with respect to the local frame of each 
fluid cell. Note that even if we use the fa, notation, fa is still given by equation (A.9 1. Also, all 3- vectors in the 
following expressions are unitary vectors. The photon's direction is given by 

p^y = (sin 9-y cos fa , sin 9 1 sin fa, , cos 9 1 ) (A. 14) 
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We need to specify the origin for (f> a (<fi a — 0). Since the (j> a integration is over 2ir, the origin is arbitrary. Let's choose 
the origin to be in the z-p 7 plane. Then, 

Pa{4>a = 0) = (sin(0 7 - 6> a )cos</> 7 ,sin(6» 7 - 9 a ) sin </> 7 , cos (6> 7 - 9 a )) 
= ( (sin 9 1 cos 9 a — cos 9 1 sin 9 a ) cos 7 , 

(sin# 7 cos# a — cos0 7 sin# a ) sin0 7 , cos0 7 cos# a + sin# 7 sin# a ) (A. 15) 

To have an expression for p" a , we need to rotate p a (<fia — 0) by an angle of <p a around p 7 . The rotation matrix is then 
given by 

" -rn r 12 r 13 

i? = ( r 21 r 22 r 23 | (A. 16) 

^31 7"32 r 33 



ru = 
ri2 = 
ri3 = 
r 2 i = 
r 22 
r 23 
r 3 i 
r 32 

r 3 3 



cos 4> a + sin 2 9 1 cos 

sin ^ 7 (1 
sin # 7 cos 9 1 cos </> 7 (1 — cos 



■ 2 </> 7 (l - COS0 Q ) 

3 # 7 sin </5> a 



So our expression for p a is finally 



sin 2 9~ ( cos 7 sin 7 (1 — cos </> a ) — cos # 7 f 

# 7 cos 9 j cos </> 7 ( 1 — cos <j) a ) + sin # 7 sin ^ 7 sin </> a 
= sin 2 # 7 cos 7 sin </> 7 (l — cos </> a ) + cos 9 1 sin </> a 
= cos (p a + sin 2 6> 7 sin 2 7 (1 — cos a ) 
= sin 9 1 cos # 7 sin § 1 ( 1 — cos <p a ) — sin # 7 cos 7 sin </> a 
= sin # 7 cos # 7 cos </> 7 ( 1 — cos (f> a ) — sin 9 1 sin </> 7 sin </> a 
= sin 9 1 cos # 7 sin § 1 (1 — cos a ) + sin # 7 cos 7 sin </> a 
= cos <p a + cos 2 # 7 (1 — COS(j) a ) 



(A.17) 



(sin 6* 7 cos 9 a — cos # 7 si 
p y a = R (sin 9 1 cos # a — cos 9 1 s: 
cos # 7 cos 9 a + sin 6 

Note that all the components of p a are proportional to either one or : 
we can rewrite, for instance, the x component as (with the A,B and C < 

p x a = cos 9 a I cos 4> a A% + sin 4> a B x + C x j + sin 9 a j < 

and in a similar way we have 



(sin 9j cos 9 a — cos 9 1 sin O ) cos </> 7 
(sin 9 7 cos # a — cos 9 1 sin a ) sin </> 7 
cos 7 cos # a + sin # 7 sin # a 

either one or no power of cos# a , sin# a , co: 
he A,B and C coefficients that only depend on photon 

(A.19) 



(A.18) 

6 a , coscj) a and sin <p a , so 
end on photon angles) 



cos 4> a A s x + sin (f> a B s x + C' x } 
C v} 



Pa = C0S 6 a 

and in a similar way 

p v a = cos(9 a { cos4> a A c y + sm<f> a B c y + C*} + sin6» a { cos0 a A; + sin0 aJ tf° + U" y j- 

= cos^ a | cos0 a A^ + sin0 aJ B^ + C^j + sin# a { cos0 a A* + sin0 a B| + } 

Now, going back to the expression for the rate, only the viscous corrected momentum distribution 
depend on <j>i so let's look at the </>i integration of the viscous corrections 

?i) + 6f(pi))(fo(E 2 ) + 6f(p 2 )) 

-fo(Ei)6f(p 2 ) + f (E 2 )5f( Pl ) 

; of 4> a , what we 



where the term p: 
need to look at is 



s / #i/o(^i)/o(^2) + 
Jo 

has been neglected. Since fo(E a ) is independent i 



(A.20) 
(A.21) 
functions actually 



proportional to Sf(pi)Sf(p 2 
is 

/•27T P 2 7T 

/ dfadfipa) a / — 
Jo Jo V 



(A.22) 
really 



(A.23) 
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It should be pointed out that 4>2 (which is given by (j>+ or 0_) differs from 4>x by a term which, as far as the integration 
over <f>i is concerned, is a constant. Because the only dependance on <f>2 comes from sine or cosine terms, and because 
the 4>i integration is over a full cycle (from to 2n), (j) 2 is actually the same as fa. Or, more precisely, 



d(f>if(cos fa, sin fa) 



dfaf (cos </> 2 ,sin fa) 



dfa f (cos 4>+ , sin fa\ 



dfaf (cos fa. , sin fa.) (A.24) 



One can replace the [j(E 2 , 2 , 4>+) + KE2, @2, fa-)) term in (A. 10 1 by 2/(i?2, #2, </>i). Furthermore, the 0i integration 

will then be treated the same way regardless of the particle being considered (1 or 2). Considering, for instance, the 
Pip\ term and carrying out the fa integration, one obtains 



2tt 



dfaPiPi 



X V 



j d(j)i ^cos^i{ cos + s\n4> a B x + C x ) + sin#i{ cosfaA s x + s\nfaB x + C x ) 

(cos6>i{ cos^iA^ + shifaBy + C y } + sin6>i{ cosfaAy + sinfaB y ' + C*} 
os 2 1 (A%A c y + B c x B c y + 2C c x C c y ) + sin 2 9, (A x A s y + B s x B s y + <2C%C\ 



costal sm^i 



(a c x a; 



A%A% 



B x B y 



(A.25) 



If we were dealing with particle 2, the #i's would be 2 's, and if we were looking at other x-y-z combinations the A,B 
and C coefficients would be different, but the general form is always the same. The important thing is that, after 
the fa integration, the p"p^ terms can easily be split into a part that depends only on E and T, the photon energy 
and the local temperature (via the cos# a and sin# a ), and one that depends only on the photon angles (the A,B and 
C coefficients). So in the end, the photon production rate, including momentum distribution function modifications 
from viscosity, is given by 



E 



d 3 R 



1 



d 3 p 16(2n) 8 E 



dt 



yJaE\ + 2bE 2 + c 



dEi 

fo(Ei)f (E 2 ) -f 



dE 2 \M\ 2 (l + f (E 1 +E 2 -E)) 
V 1 



s 2T 3 



[/o(£ 2 )/o(£i)(l + /o(£i)K^ 



f (E 1 )f (E 2 )(l + f Q (E 2 ))p%p% 



(A.26) 



where it is understood that the p^Pa terms have been integrated over fa, from to 27r. Now if we substitute Eq. 
(A.25 1 (and similar expressions for other p^Pa combinations) into the above expression, the result can be expressed 
in terms of four-integrals for which the result depends only on E and T, multiplied by a combination of terms that 



depend only on 9 1 and 



ds 



where [...] is one of 



The integrals that need to be computed are then of the type 



dt 



dE x 



dE 2 \M\ 2 (l + f (E 1 + E 2 -E)) 



1 



Va~Ef 



2bE 7 



-.[...] 



(A.27) 



[...] = 



' h(Ei){E 2 2 


-ml) 


cos 2 


02/o(£ 2 )(i + /o(£ 2 )) 


h{Ei){E% 


- rn\) 


sin 2 


^/o(S 2 )(l + / (^2)) 


fo{Ei){El 


-ml) 


cos 6 


2 sm0 2 / o (£ 2 )(l + /o(# 2 )) 


< 

HE 2 ){El 


— m\) 


cos 2 


6 1 f (E 1 )(l + f (E 1 )) 


ME 2 )(E 2 1 


-m\) 


• 2 
sin 


9ifo(Ei){l + f (Ei)) 


> ME2KE 2 


- m\) 


cos 6 


1 sin0 1 /o(^i)(l + /o(^i)) 



(A.28) 
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where (E% — m 2 ) is just \p a \ 2 , which we have omitted in the derivations above. For instance, to co mpute the 
contribution to the correction from the pfp\ term, one would need the last three terms in Eq. (A.28) (the ones 
corresponding to cos 2 6* x , sin 2 0\ , cos 0% sin 0\). Each of these terms, would then be multiplied by some combin ation o f 

the A, B, and C terms. In this example, the cos 2 6\ gets multiplied by [A c x A c y + B c x B c y + 2C x C^j , as in Eq. ( A.25 ) 
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